o 
o 



0^ 
O 



Energy-Based Ghost Force Removing Techniques for the 
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Abstract 



This paper studies numerical methods for accurate treatment of the interface between the 
local and the nonlocal region in a QC approximation of atomistic materials. Only the energy- 
' based methods are considered. Particularly, a quasicontinuum projection (QCP) method based 

I on the idea of finite elements is shown to be accurate and efficient for this problem. We analyse 

the QCP method and study its relation to the existing methods, such as the quasinonlocal quasi- 
continuum method and the geometrically consistent reconstruction-based method. The analysis 
^ • and the results of numerical tests confirm that the projection-based QC method successfully 

removes the ghost force with the same computational cost as the other methods. In all com- 
puted examples the error of QCP is either the same or lower as the error of the other methods. 
The performance of these methods in treating interfaces of elements in the local region is also 
^ ' examined. 
I> 

m 

^ ■ 1 Introduction 
0\ _ 

I Many fundamental problems in science and engineering can be modeled by partial differential 

equations, which rely on the macroscopic or continuum description of the problem. In spite of 
tremendous success of continuum models, they also have certain limitations. For example, when 
modeling macroscopic materials with microscopic defects, a continuum model will be less accurate 
^ ' near the defect and may give even wrong predictions, because the continuum model neglects the 

■ microscopic features of the defect. Therefore, one might be tempted to switch to a complete 

atomistic model that takes the full atomistic description into account. However, this is not an 
efficient strategy of modeling macroscopic materials not only because large atomistic systems are 
too large to handle even with the existing most powerful computers, but also because the solutions 
we obtain will likely contain too much information that is of little interest. Thus, a combination 
of both atomistic and continuum descriptions is required for successful modeling of materials with 
defects. Such an approach is called atomistic/continuum multiscale modeling. 

Atomistic/continuum models adopt the point of view that there is an underlying atomistic 
model of the material which is the "correct" description or exact solution of the material problem. 
This could be a quantum-mechanically based description such as density functional theory, but in 
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practice it has been based primarily on classical-mechanics atomistic models with semi-empirical 
interatomic potentials. The problem thus has the two length scales: the length of the whole material 
and the typical interatomic distance. The ratio between these scales, e, is usually much less than 
1. The bulk behavior of the material (at least in the case of small deformations) can be described 
by PDEs of continuum mechanics. However, if one needs to model the material defects of the size 
0(e), then one may require employing a fully atomistic model near those defects. 

A representative of atomistic/continuum models is the quasi-continuum approximation (QC) 
which is becoming a popular multiscale technique for simulating static properties of polycrystalline 
materials. It was put forward by Tadmor, Ortiz and Phillips [29] and later was updated in a review 
paper by Miller and Tadmor [19] . The idea of this approach is that we consider at the macroscopic 
scale the region in the material (called the local approximation region) where no serious defects 
occur and the theory of continuum material elasticity (or a coarser mesh approximation) may apply. 
In the (nonlocal) region where serious defects occur, much finer mesh or atomistic model has to be 
used. The QC method is usually described in a finite element framework. 

While a significant body of knowledge about QC and related models and their numerical suc- 
cesses has been accumulated, not much has been reported on the analysis of these models. The 
convergence rate of the QC method has been tested in [H] via computational results of a specific 
nano-indentation problem. Other numerical demonstrations can also be found in the literature, e.g. 
[271 129j . Rigorous numerical convergence analysis of the QC method has been published in [TS] for 
a one-dimensional case in the absence of external forces. Later studies with certain external forces 
have been conducted in [U [101 [HI [22] for one-dimensional and occasionally two-dimensional cases. 
An analysis of the method for two-dimensional cases under a few physically reasonable assumptions 
for solid materials has been given in [1^. Other recent developments on the QC method and other 
atomistic/continuum models from the mathematical point of view include analysis in [11] concern- 
ing the Cauchy-Born rule and in the context of the heterogeneous multiscale method [8], in [3] for 
excluding a spurious finite element effect in a prototypical multiscale model, in [9] for a necessary 
and sufficient condition of the uniform first order accuracy, in [H [21] for desirable a posteriori error 
estimators associated with the QC method, in p.8j for analysis of the cluster-summation rule, and 
in [2] for blending techniques in the mixed atomistic and continuum domain. 

A problem having being associated with the QC method since its invention is the error at 
interfaces (edges or face) of elements and at the interface of local and nonlocal regions. In [T6] the 
error at edges of elements has been estimated in local regions for QC taking representative atoms at 
the center of finite elements. Such an error increases with the number of edges, so at some critical 
point a mesh refinement would increase the numerical error even in the local QC case. The error at 
the interfaces of local/nonlocal regions may be more difficult to study since it is usually associated 
with material defects. Effort has been focused on unphysical forces (so-called ghost forces) at the 
local/nonlocal interfaces. Such forces exist in the original QC method, unless only nearest neighbor 
interactions are taken at the local/nonlocal interfaces. A correction force has been introduced in 
[27) and the method was called "force-based correction". Convergence of an iterative procedure 
based on this version of QC has been shown in [5J. Despite its ability to eliminate the ghost force, 
the method is not associated with the "correct" potential energy and thus leads to an issue with 
energy conservation and possibly large numerical error (see e.g. [28]). This is why we focus on 
the energy-based approach in this paper. The reader can refer to [Tj [20] for recent analysis of the 
force-based approach. A seamless coupling in the interface region has been sought in [28] based on 
the changing the interaction energy of certain atoms and the resulting method is known as quasi- 
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nonlocal method (QNL). A convergence analysis for the linearized one-dimensional periodic case 
has been given in ^ with the assumption of second nearest neighbor interaction. In general the 
QNL method does remove the ghost force but is limited to potentials with a relatively short range. 
In [9l [12] it was shown that QNL works up to the second nearest neighbor interactions. In pjj 
the accuracy of a more sophisticated approach which is based on what was called a geometrically 
consistent reconstruction (GCR) introduced in [9] was studied in detail. It was also claimed there 
that the GCR-based method includes QNL as a special case. All these methods consider only 
reducing errors at the local-nonlocal interface. However, the problem of reducing the error at edges 
of elements in the local region has not been addressed. 

In this paper we aim to reexamine these existing methods, particularly focusing on existing 
versions of the QC method which is not originally proposed to remove the ghost force. The method 
is based purely a finite element idea and often appears in literature (e.g. [T71 |22l [23l HI]). It is 
also briefly mentioned in [25] for a linearized model that the method does not have the ghost 
force. We shall call it the quasicontinuum projection method (QCP) in this paper. Under a one- 
dimensional setting of the problem we will carefully study relationships between QNL, GCR and 
QCP, with a particular attention on the way these methods remove the ghost force. We will show 
that QCP is not equivalent to QNL, not equivalent to the particular instance of GCR presented in 
[9], but there exist particular instances of GCR which are equivalent to QCP. In addition, the QCP 
method offers a natural way of accurate treatment of edges of elements in the local region and thus 
eliminates the error at element interfaces. This makes it an essentially more accurate method than 
the other ghost-force removing methods. Based on that we conclude that QCP is the easiest and 
the most efficient way to remove the ghost force and the element edge errors in QC. Furthermore, 
in favorable cases it does not introduce more computational costs in comparison with QNL and 
GCR. In unfavorable cases, for instance when a 2D or 3D atomistic lattice is not aligned with the 
QC triangulation, QCP can be modified in the local region (while preserving the same accuracy in 
the local-nonlocal interface) so that the number of operation is of the same order as QNL or GCR. 

The paper is organized as follows. The problem of an equilibrium of an atomistic material is 
introduced in section [2] under the two settings: the general setting and the one-dimensional periodic 
setting. It will be shown that QCP has no ghost force under the general setting in subsection 14.51 
QC, QNL and GCR, their relations to each other, and their treatment of the ghost force are 
examined under the periodic setting in section [3l QCP is then introduced and compared with QNL 
and GCR in subsections I4.1H4.4[ Numerical experiments under the both one-dimensional periodic 
setting and general setting are given in section [5] to show the performance of these methods. The 
results are discussed and summarized in section [6l 

2 Problem Formulation 

We first describe the general problem formulation and then will consider the ID period setting (as 
considered, for instance, in [6]). 

2.1 General Setting 

Consider a problem of finding an equilibrium configuration of atoms in an atomistic material in 
space Mf^. Denote the degrees of freedom (which are chosen to be the coordinates of atoms) by a 
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vector 

u = {u.i}, {l<i<N). 
To model the equilibrium of atoms, we introduce the potential energy of the atomistic system 

n(u) = ^totH + ^cxtN, (1) 

which is a sum of the internal energy of atoms Etotiu) and the external potential Ecxtiu), the later 
will be considered to be a linear functional of u (which means that the external forces on atoms do 
not depend on their positions). 

In the above notations the problem of finding the equilibrium configuration of atoms is written 

as 

^ = (i = l,2,...,iV). (2) 

OUi 

Sometimes, instead of simply an equilibrium configuration, one specifically requires a stable equi- 
librium configuration, which is formally written as 



u is a local minimizer of Il{u). 



(3) 



The equilibrium equations ([2]) are often solved by Newton's iterative method which in matrix 
form reads 



u 



duiduj 



u 



n+l 



U 



an('u") 

dui 



To assemble this system of equations one should effectively compute the right-hand side vector 
components and the stiffness matrix components . The algorithm of computing these 

would depend on a particular problem setting. Below we elaborate the ID periodic setting, since 
we find it to be the easiest formulation for analysis and comparison of different QC approaches, 
especially if one's focus is the ghost force. 



2.2 ID Periodic Setting 

In this subsection we describe the problem formulation of finding an equilibrium of an atomistic 
material in the ID periodic setting. We consider the periodic boundary conditions, as it is done in 
[6] . This is done to avoid difficulties arising from presence of the boundary of the atomistic material. 
Otherwise, the boundary of an atomistic material, unless properly treated, would contribute an 
additional error to the numerical solution. 

Consider an atomistic material in one dimension with atom positions Uj, (— cxd < i < oo). We 
assume that the material behaves periodically with a period of length one over atoms: 

Ui+N = Ui + 1 (— OO < i < oo). 

Thus, the degrees of freedom of such system are described by 

u = {ui}, il<i<N). 

The energy of atomistic interaction of the system (summed for the atoms over one period) is then 

i = l j = — OD 
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where ip{r/e) is the potential of interaction of atoms at distance r. For compactness of notations, 
we define the potential if{z) for negative z as ip{z) = ip(—z). We assume that the potential if{z) 
vanishes for \z\ large enough, so that we need to consider at most n neighboring atoms in the 
interaction energy: 



N i+n , ^ 
%=\ j=i—n ^ ^ 

The potential energy of the external force / is 

N 

Ee^tiu) = -e'^fiUi. (5) 



i=l 

The forces fi on each atom are given and considered to be independent of atom positions Ui. For 
the problem to be well-posed, the sum of all forces is assumed to be zero: 

N 
i=l 

If the sum is not zero, then the total energy will not be bounded from below due to the periodicity. 



3 QC approximation 

Usually, the characteristic distances on which the external force / varies are much larger than the 
characteristic interatomic distance e. In this case we do not need to model all degrees of freedom of 
the system associated with each individual atom and can employ a QC approximation. It consists 
in choosing the nodal atoms Inod = {h,'i'2, ■ ■ ■ ^'^k} and regarding the positions of the non-nodal 
atoms to be known via a linear interpolation: 

Ui = .^^ .'^ Ui^_^ + ^ — ^-^^^Ui^ {ik-i < i < ik), k = l,...,K. (6) 
— ik-l ^k — ifc-l 

Here we defined io = ^i^ ~ by periodicity. Then all distances between neighboring atoms posi- 
tioned between ik-i and ik are equal to the same value 

Ui - Ui-i = — : [Ik-l < I < Ik)- 

^k — l-k-l 

Thus, the unknowns in a QC approximation are 

u^'^ = = [uij^. (7) 

After having defined the nodal atoms, we need to approximate the original problem The 
original QC method (see e.g. [HI [19l [29] ) is based on the assumption that the energy of the segment 
ik-i ^ i ^ ik can be approximated by the continuum energy via the so-called Cauchy-Born rule: 
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^" • nodal atom 



/; ii i2 ... iK Ii ' non-nodal atom 

Figure 1: Illustration of local and nonlocal regions 



periodic extensiony v Periodic extension 

' • I • • I : • • 

'0 ^ 'I '2 ••• IK ^ 'K+\ 

Figure 2: Illustration of periodic extension 



where the interatomic energy density is defined as 

^ n n 

^(^) = 2 M = ^ V9 (mz) . (8) 

m=—n 771=1 

In the present paper, we are focused on modeling a localized irregularity of the atomistic mate- 
rial, which requires keeping track of individual atoms in the neighborhood of the irregularity (called 
nonlocal QC) while using a local QC approximation elsewhere. To model such situation in ID, we 
consider the external force / to be zero everywhere except for the atoms M, M + 1,...,M + P — 1: 

/i = (i = l,...,M-l,M + P,...,iV), 

where M is chosen to be M = [A^/2j. 

A common approach to this problem is the following [26]. Introduce the nonlocal region 

T^ = {i: M -L<i<M + P-l + L} 

and the local region 

Ii = {l,2,...,N}\I^. 

The nonlocal region is chosen as an extension of the region where fi^OhyL atoms in each 
direction. The nodal atoms Inod = {h,i2, ■ ■ ■ lix} are chosen so that Xn C Inod- Since we assumed 
the external force / to be equal to zero in the local region Ii, the deformations will be well-defined 
by a linear function and therefore nodal atoms are not required to be added inside 2i (see figs. [Hand 
[2]). We stress again that this specific problem setting is chosen to study the relationships between 
ghost-force removing strategies. In a more general setting, local and nonlocal regions, as well as 
the mesh refinement in local regions, may be determined through a certain adaptive strategy (see 
e.g. [El [29]). 

Thus, the expression of the total internal energy of the standard QC method is 

T^QCE T^QCE . T^QCE 

-^tot — -'^local "T -^nonlocal' 

where 

\e{ii-io)J 
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for $ defined by ([8]), and 

p^QCE ^ 
' ' ' 2 



iK i+n 

QCE iy y ( ^^-^3 

nonlocal 9 / j / j " 



i=i\ j=t—n 



The respective system of equations to be solved is 

anQCE ^ ^^QCE gj^^^^ 



+ir^ = (fe = l,2,...,K), (9) 
dui^ dui^ dui^ 

which is the shortened version of the original system ([2]). 

3.1 Ghost Force, QNL Method, and GCR method 

The traditional QC method for atomistic material modeling, though seems to be very intuitive, has 
an essential drawback: the accuracy of approximation in the neighborhood of the interface between 
Tn and Z\ falls down to zeroth order. This phenomenon of loss of accuracy on the interface between 
the regions is commonly interpreted as the ghost force (see e.g. [H 119^ [26]). 

The nature of the ghost force can be illustrated in the case of second nearest neighbor interaction 
n = 2. Consider the equidistant lattice Ui = ezi with the zero external force / = (and hence zero 
E'ext)- In this case the force on each atom should be zero (due to reflection symmetry). Compute 
the force on the atom 12'. 

«rQCE jifQCE fiirQCE o «K i+n / \ 

dE^ot _ ^-^nonlocal , ^^^local ^ ^ ^ ( ' ^J ] , n 

f),,. Of),,- ^\ ^ 



dui2 dui2 dui2 2 dut.2 



z=ii j=i—n 



which upon substitution Ui = ezi yields 

The reason for not getting zero force is that we are missing the term 

|(^((Mi2_2 -UiJ/e) 

in the expression for E^^^ (specifically in E'^^). 

There are several methods to eliminate the ghost force, most of which can be divided into the 
two categories: force-based methods (which modify the right-hand side of the equilibrium equations 
dl])) and energy-based methods, which modify the interaction potential e'^^ . As mentioned in 
Introduction, we will not discuss any force-based methods but focus only on the energy-based 
methods. One of the energy-based methods proposed to eliminate the ghost force was the so-called 
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quasinonlocal quasicontinuum method, or QNL in short [28]. The idea of the QNL method is to 
modify the interaction of atom 12 (and the respective atom on the other side of the nonlocal region 

iK-i) in order to eliminate the ghost force. Namely, instead of the term ip 



the term ip 



2{ui2-Ui2-i) 



we introduce 



I , just as in the Cauchy-Born extrapolation, but for the nonlocal atom ^2. 
In this case the atom 12 is called a quasinonlocal atom, since though it formally lies in the nonlocal 
region, it interacts with atom Z2 — 2 like a local atom. Thus, the expression of the energy of the 
QNL method is 



pQNL pQNL , pQNL 



tot local ' nonlocal ' quasinonlocal' 



where 



for $ defined by ([8]), and 



Ui 



Ho 



e{ii - iq) 



e: 



■QNL 
nonlocal 



■ E E 

i=ii j=i—n 



Uj — U-, 



and finally 



^QNL ^ £ 

quasinonlocal 2 



2(u 



12 



^22 — 1 



1ij2 ^42 + 1 



u, 



«2 



-2>^ 



iK-l 



Ui 



The QNL method has been designed to have no ghost force for the case n = 2, which means 
that the force exerted on the atoms in the equispaced lattice Uj = ezi is exactly zero. One, however, 
can do a more detailed analysis: consider an arbitrary lattice satisfying the QC constraints ([6]) and 
compare the QNL energy with the exact atomistic energy: 



-^tot - ^tot - 2^ 



+ 2^^ 



+ respective terms for atoms around ix- 

Defining D'^Ui = 2m;+Hi_i expanding E^^ — Etot in Taylor series w.r.t. Ui^ around its 
extrapolated value 2ui-^ — Ui-^-i (the respective terms for atoms around ix are expanded w.r.t. 
Uij^-i) yields 



pQNL p 



+eO — ^ + respective terms for atoms around 



IK 
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Ui 



+ eO [eD^Ui 



) 



Ui 



+ eO ieD^Ui 



One can see that if the interface atom ii (as weh as ik) hes in the region where the solution is 
smooth (i.e. D'^Ui^ is bounded, then the difference E^^ — Etot is small. Moreover, the force on 
the atom 12 is now 



dE: 



iQNL 
tot 



dEi 



tot 



dui 



duj. 



'-{D\,,f + 0{e^D\.,,f), 



"12 •^"■12 

which is 0(e) if the solution is smooth near the interface. In the next section we will show that 
the QCP method has exactly the energy -Etotj which means that it is not equivalent to the QNL 
method. 

E, Lu, and Yang [9] proposed a way to remove the ghost force for an arbitrary interaction 
distance n, based on what they call a geometrically consistent reconstruction (GCR). They inter- 
preted the change of the interaction of certain atoms as changing the reconstruction of atoms: for 
instance, in the QNL method, in the interaction of atom 12 with 12 — 2, the position of atom i2 — 2 is 
reconstructed as + 2{ui^ —^12)- The idea in [9] was to change the interaction (i.e. reconstruction) 
of more than one atom near the interface (the number of atoms whose reconstruction is changed 
depends on n) so that the ghost force does not appear for a given n. The reconstruction was 
sought by the indeterminate coefficients method in the form of a linear combination of the nonlocal 
reconstruction and the Cauchy-Born reconstruction: 



u 



reconstructed 

'i 



c: 



GCR 

ij 



Ui + (1 



Cij*^^) {ui + (i - ^)(^i^+sgn(i-^) - Ui)). 



(10) 



The requirement of such reconstruction to have no ghost force results in a class of methods. E, Lu, 
and Yang gave particular instances of such methods in [9j . For n < 3 the coefficients Cij given in 
table I in P p. 214115-8] are: 



GCR 



' 1 


(iJ) 


= {ii- 


l,ii) or 


{iK + l,iK), 




1 


(iJ) 


= in- 


l,ii + l 


) or {iK + l,iK 


-1) 


1 


ihj) 


= in- 


l,ii + 2 


) or (iK + l,iK 


-2) 


2/3 


ihj) 


= i^l- 


2,ii + l 


) or (iK + 2, ix 


-1) 


1/3 


(iJ) 


= (n + 


l,k-2 


or {ix - l,iK 


+ 2) 


/^QCE 


otherwise. 









where the coefficients C^^^ are defined as follows: 



a 



QCE 



1 i < ii or i > iK, 

ii < i < iK, 

1 « € {n, iii"} and ii < j < ix, 

« € iii"} and {i < ii or i > i^) 



For example, as is implied from table I in [9, p. 214115-8], for n 
brevity written below as a correction to E^^^) is 



TP' 



GCR 
tot 



T-iQCE . 
^tot + 



-1 — "Uii+l 



2(ni,. 



2 the atomistic energy (for 

-Ui,)^ 
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e 




Note that here we ignored another correction term 




e 



) 



- 9? 



( 



e 



which equals zero identically under the QC constraints but is implied in the table in [9j. It is 
worth noting that this particular instance of the GCR method is not equivalent to the QNL method 
for n = 2 (although [9] claimed it is equivalent): the QNL does not change the interaction of the 
local atom ii — 1, whereas that particular GCR method does. However, it is true that the QNL 
method is a particular instance of the GCR method. Particularly, if the one changes the coefficients 
to C^-^ in table I of [9j for n = 1 and n = 2, then one will obtain exactly the QNL method. 
That is, in notations of this paper the QNL method's reconstruction is 



which is clearly different from C^j . 

Another property of the instance of the GCR method in table I of [9] which is worth noting 
is that atomistic energy of the GCR method identically equals to the exact atomistic energy: 
-^to?^ = ii^tot for n = 2 for a lattice satisfying 1^. As we will see below, it is not the case for n = 3, 
although one can find other instances of the GCR method, which would give the exact energy E'tot 
on the QC lattice for a given n. 

GCR gives a general framework or criterion for removing the ghost force at the local-nonlocal 
interface. Nevertheless, it involves a priori calculation of reconstruction parameters and the for- 
mulation of the method is different for any slight change of the problem, e.g. from n = 2, 3 to any 
given number. In the present paper we thus seek and study a specific method simply based on the 
projection ([6]) requiring no a priori calculation of any parameters and meanwhile performing not 
worse than GCR in all occasions. We will call it the QC projection method (QCP). 



The idea of projection was used by Rudd and Broughton in [23|, \2^ [25] in their coarse-grained 
molecular dynamics method and the absence of the ghost force for such method was noticed by the 
same authors in [25] for a linearized model. A rather detailed analysis of such method can be found 
in [22] for a one-dimensional steady-state problem satisfying the Dirichlet boundary condition. In 
this method, the full atomistic solution is "projected" onto the space of admissible deformations 
([6]), in the same sense as the solution to differential equations is projected onto the finite element 
space in Galerkin projection methods. We therefore call it the quasicontinuum projection method 
(QCP). This approach can also be interpreted as introducing the constraints ([6]) to the original 
optimization problem ([2]). In this case the approximate solution u^'~' (cf. ([7])) will be the best 
possible QC solution, i.e. the solution which minimizes the energy n(it) among all u satisfying ([6]). 

In this section we first formulate the quasicontinuum projection (QCP) method in the general 
setting, then we formulate it in the ID periodic setting for the convenience of discussing the 
relationship of QCP with the other ghost-force removing methods near the local-nonlocal interface. 




otherwise 




4 QCP method 
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4.1 QCP in General Setting 

We want to approximately find the critical point of !!(«) by reducing the number of degrees of 
freedom of u by considering the constraints ([6j). A natural approach would be to just restrict our 
search for the critical point to the subspace defined by the constraints 
Formally this can be written as follows. Define the QC variables 



u 



QC 



which are analogous to those in ID periodic setting ([7]). Define the QC reconstruction matrix 
in the following way: 

u = C/QC^QC^ (11) 

which corresponds to the QC reconstruction ([6]). Then the QCP method can formally be written 
as 



d 



dui 



n {u^'^u 



{l<i<K), 



(12) 



^k 



or alternatively, in relation to the stable equilibrium problem ([3]) as 

u^'~^ is a local minimizer of IT {U^^u^^^. 
Define the QCP approximation to the potential energy as 

nQCP (^QC) ^ n (t/Qc^QC) . 

The equilibrium QC equations ([12]) can also be solved using the Newton's method 



{u 



QC,n+l _ ^QC,n\ 



kl 



duk 



(13) 



(14) 



The implementation of the QCP method in the general setting will be discussed after further 
elaboration of the QCP method in the ID periodic setting. 



4.2 QCP in ID periodic setting 

In this subsection we formulate the QCP method in ID periodic setting and analyze it. For that 
we take the specific form of the reconstruction operator [/Q'-' which corresponds to formula ([6]): 



u 



QC^ 



Ik - I I- lk-1 
: + : U 



Ik — ^k-l 



ik — ik-1 



Here we continue using the periodic extension of indices: ix+i = h and io = ix- Notice that the 
matrix [/Q^ is sparse: each row of [/Q"-^ consists of at most two entries. 
Now, according to p2|) . we need to write down the potential energy 



u 



QC^ 



n([/Qc 



u 



QC^ 



which is nothing but the energy ([T]) written in terms of positions of the nodal atoms Uj^ , . . . , . 
The external energy E'ext is easy to deal with since / is zero for the non-nodal atoms. To write 
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down the atomistic interaction energy of QCP, we split it into the three parts, corresponding to 
interaction between nonlocal atoms {ii < i < ir), interaction between local atoms (i < ii or 
i > ix), and interaction of nonlocal atoms with local atoms: 

^QCP ^ ^QCP ^ ^QCP ^ ^QCP _ (^5) 

The energy of interaction between nonlocal atoms E^^^ can easily be written down since all 
nonlocal atoms are nodal atoms: 

0<\i-j\<n 

The energy of interaction between local atoms E^'^^ must be written in terms of the nodal atoms 
and Uij^. For this purpose, we denote iq = ix — N and will use the periodicity condition 
Uj-jv = Ui — 1 for writing down this energy. We will also use Ui — Uj = (^ — ■ Thus, the 

energy E^ is 

fQCP ^ £ V- / -Uj^) 
2 2 . . ^\ e{H-io) 

0<\i-j\<n 

m=—n ^ ^ / (17) 



n 

= e ^ (ii - io + 1 — ''T^)<^ 

m=l 



e(n - io) 

(«l - «o)¥' 7^ - e >^ ("^ - 1)<^ — T. 



m=l ^ m=l 



Here we see that the first term in the final expression for E^^^ is exactly the expression of the 
Cauchy-Born rule: Ylm=i f {' ) is the strain energy density associated with a representa- 

tive atom and e(ii — io) is the length of the strained material. We can also see that in the case 
of nearest neighbor interaction (n = 1) the Cauchy-Born rule gives exactly the local component of 
the QCP energy and hence no corrections are required. However, if n > 1 then one needs to add 
the third component of energy (in order to avoid the ghost force), namely the energy of interaction 
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of local and nonlocal atoms (atoms in the neighborhoods of ii and ix)- 



e: 



QCP 



j<h<i 
0<|i-j|<n 



i<iK<i 
0<|i— j|<n 



Ui — Ui 



j<ii<i 
0<\i-j\<n 



(18) 



Ei'ui - Ui^ (ii - j){ui^ - Ui^) 
^\ + 



e(ii - io) 



i<iK<j 

o<l«^il<" 



+ 



e(«x+i - ix) 



Uu + 1- 



Here we denoted ir+i = ii + N and implied u 

Thus, we have written the interaction energy in terms of our unknowns . . . ,Uij^. Adding 
it up with the external potential energy will yield the total potential energy of the QCP method 
HQCP rpj^g solution is defined critical point II'^'-"^. Finding a critical point typically 
consists in starting with some initial guess and performing Newton's iterations ()14p . For Newton's 



iterations one needs the right-hand side vector 



and the stiffness matrix 



, both 



kl 



of which can be computed from the above formulae in a rather straightforward manner. We discuss 
the implementation of the projection method in subsection 14.41 

First, note that in the case of nearest neighbor interaction (n = 1), E^'^^ is exactly the 
continuum energy computed by the Cauchy-Born rule and E^ = 0. It means that the QCP 
method is equivalent to the traditional QC method in this case. 

OOP 

Second, note that the local component of the energy E2 does not include interaction of 
nonlocal atoms with certain local atoms close to the interface (namely, local atoms ii — n+l, . . . , ii — 
1 and i/^ + 1, . . . , iii: + n — 1). These atoms play the role similar to the quasinonlocal atoms in the 
QNL method [28] , with the difference that in the original QNL method only nonlocal atoms change 
the way they interact, whereas in the QCP method the local atoms close to the interface change 
the way they interact. The QNL method successfully removes the ghost force for second nearest 
neighbor interaction {n = 2), but still exhibits the ghost force for longer interaction distances 
(though the magnitude of the ghost force of the QNL method is less than that of the QCE method, 
as we will see from the numerical results presented in section ^ . 



4.3 Comparison of QCP with GCR 

In this subsection we will show that, in relation to the local-nonlocal interface, QCP is not equivalent 
to the particular instance of GCR that was presented in [9], but one can find an instance of GCR 
which would be equivalent to QCP. However, GCR and QCP are essentially different in the local 
region in the way they treat interfaces of elements in the local region: GCR (as well as QNL) 
uses the Cauchy-Born extrapolation as in the standard QC, while QCP uses the projection. The 
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Cauchy-Born extrapolation and the projection give different results in the case of a non-equidistant 
lattice in the local region, since they treat interface between elements in a different way. As was 
shown in [16], the Cauchy-Born extrapolation introduces an additional error at element interfaces 
in the local region. This additional error may be crucial, as we will see in subsection 15.21 

One characteristic feature of the QCP method is that the QCP atomistic interaction energy 
coincides with the exact atomistic energy for the atomistic lattice satisfying the QC constraints ([6]). 
It therefore is not equivalent to the QNL method for n = 2 (see the calculation in section 3.1 which 
shows the QNL energy is different from the exact atomistic energy) and therefore not equivalent 
to QNL for any n >2, but is equivalent to GCR for n = 2 as mentioned at the end of section 3. 

For the case of n = 3, substitution of the values of coefficients for GCR from table I of [9j and 
direct computation (which is straightforward, but too bulky to be presented here) yields 

pGCR pQCP e f 2{ui^+i - Ui^^2) , 



3e 3e 



e f 2{ui^_2 - Uii+i) , Uii-2-Uii-i 

^ \- 



2' \ 3e 3e 

-e(p ( ^lil ij_2 j ^ respective terms for atoms around ix- 



As before, defining D^ut = a,nd expanding E^f^ — E^^^ in Taylor series w.r.t. ut^ 

around its extrapolated value 2ui^ — Uj^-i (the respective terms for atoms around ik are expanded 
w.r.t. Uij^_^) yields 

One can check that the atomic force —g^ — equals to up to the first order of accuracy in e 



provided that the underlying solution is smooth (i.e. that D m is bounded). For instance 

Moreover, the quantity rather small in practice: for instance in the case of 

Lennard-Jones potential, if the interatomic distance is approximately equal to e, then 

2 / 3(n,-. ) \ ^ 2 ^ 
9^ V e{H-io) J 9^ ^ ' 

Therefore we should expect that the difference in solutions by the QCP method and the GCR 
method is small, provided that the solution is smooth around the interface between the nonlocal 
and the local region (for D^Ui not to be too large). However, it should be noted here that if the 
solution is not smooth enough in the local region, then the QCP method may give significantly 
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better results as compared to the GCR method (since they treat interfaces of finite elements in the 
local region in a different way). This will be shown in the section with numerical tests. 

Having established that the instance of the GCR method presented in table I of [9j is not 
equivalent to the QCP method, it is nevertheless possible to come up with the GCR method which 
would be exactly equivalent to the QCP method. For that, we should obtain the GCR method 
with the nonlocal atoms reconstructed as fully nonlocal atoms. That is, the corresponding C^^^ 
in (jlOp should all be equal to one if i is the nonlocal atom, except for the interface atom (ii or ix 
in our case) which is allowed to be less than one if j is the local atom. This can easily be attained 
by shifting the rows of table I to right by n ^ ^ o „ u,, 

-(n-l) 



1 for n > 2 (i.e. by assigning new values of to the 
old values of C^_^^^_^-^ in notations of [9]; note that for n <2 the method is already equivalent to 
QCP). 

Computations below verify that for n 
coefficients C^^^ as described above) 



3 the following GCR method (obtained by shifting 



a 



GCR 



1 
1 
1 
1 

2/3 
1/3 

(^QCE 



(i,i) = {ii - l,ii) or {iK + I, Ik), 
{hi) = {h - l,n + 1) or {iK + l,iK - 1), 
{hi) = {h - 2,ii + 1) or {iK + 2,iK - 1), 
= in -l,k + 2) or {iK + l,iK - 2), 
= {ii - 3,ii) or {ix + 3,ix), 
(ij) = ih,k - 3) or {iK,iK + 3), 
otherwise 



(20) 



is equivalent to QCP. Indeed, the internal energy for the QCP method for n = 3 is 



pQCP _ pQCE 
-^tot — -'^tot 



+ 



2^^ 





-1 




Ui^ + l 






e 






-2 




Uil + l 






e 






~1 







2^ 



e /2(uj^_i-niJ 

2^ ^ 

e /3(-Uj^_2 - Mii_i) 
2^ 



3(Uj^_l - Mi J 



e J z \ e 
+ respective corrections for atoms around ix- 
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The formula for the energy of this GCR method, according to (j20p is then 



pGCR _ pQCE 



+ 
+ 
+ 



+ 3 2^ 



3 V2 





e 


■"11-2 


- -"n+i 




e 




- '"n+2 




e 








e 


( 


- Uh-3 



2^ 



2^ 



2{ui^ 






-1 - ^ijy 


3(^11 


-2 - ^iii-l) 


3(^11 


e 

-1 - 


3(-Ujj 


-3 - 




e 


3(iiji 




e y y 



(21) 



+ respective corrections for atoms around ix- 



Notice that the term in parenthesis in the first hue of ()2ip is identicahy equals zero (this term 
was kept to match the formula (j20p ). Taking into account that these energies are computed on an 
admissible QC lattice satisfying the constraints ([6j), one can also notice that the terms in the fourth 
and the fifth line of (|2ip are identically zero, which makes E^^ = E^^ for n = 3. However, as 
we mentioned earlier, QCP and GCR are different if there are nodal atoms in the interior of the 
local region. The subsection 15.21 shows that this may be crucial in some cases. 

We would like to summarize the analysis of the QCP method with one important remark. In 
estimates such as (I19p . following the conventions set in earlier works, we estimated the QC energy 
or the QC force through the non-nodal atoms such as One should understand this in the 

following way. The QC energy, as introduced in this work, depends only on positions of the nodal 
atoms u^'^ . Therefore, strictly speaking, the positions of the non- nodal atoms should be expressed 
through u^'^ using ([6]) and only then used in the estimates. For instance, one should express Ui^-i 
as 

- ^ , h-i^- 1 

n - iQ ii- io 

Thus, the estimate ()19p in its expanded form should read 



iP 



3(nj^ 



e{ii - io) 



where D'^Ui^ should be expressed as 



H2 



Ui 



4.4 On Implementation of the Projection Method 

This subsection is devoted to the implementation of the projection method. We first describe 
the possible implementation in the ID periodic setting and then discuss the implementation in 
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the general case. The implementation of the projection method consists in computing the right- 



hand side vector 



92nQCP 



, which are the first and the second 

kl 



and the stiffness matrix 

k 

derivatives of the total energy with respect to the positions of the nodal atoms Ui^, . The total energy 
of the projection method is n^^^ = E^^^ + E^^ , where E^'^^ is defined by equations (fT5|) - p^ 
and E^^ = E^xt is given by ([5]). Computing derivatives of E^^ is fairly straightforward since 
it is a linear functional of our unknowns Uj^,. That would be somewhat less straightforward if we 
assumed a smooth nonzero external force in the local region: in this case we need to express -Eext 
through the positions of the nodal atoms Uj^. and possibly use the approximate nodal summation 
techniques to optimize the performance of the method. 

Differentiation of atomistic energy E^ is relatively straightforward: one should go through 
all atoms that contribute a nonzero interaction potential with the given atom and compute the 
respective entries of the matrix and the right-hand side vector. Computation of derivatives of 
E2 is also staightforward: it is the standard QC energy with some correction terms. The energy 
S^^^ contributes only to the entries associated with the nodal atoms in the local region {ii and 

OCP OCP 

iK in our case). Computation of derivatives of E'g is somewhat less straightforward since -Eg 
involves the interaction energy of the nodal and the non-nodal atoms. 

OOP 

The strategy of computing derivatives of E^ is the following: we go through all pairs of 
atoms Ua, Ufs that give a nonzero contribution to the interaction energy, where the atom Ua = Ui^. 
is in the nonlocal region and the atom up is in the local region. For the atom in the local region, 
we should represent its coordinate through our unknowns nj^: 

up = Ul [ui^\ , 

where [u^]l is the vector of unknowns and Up is the vector of coefficients of representation of up. 
Note that the vectors Up are essentially the columns of the matrix C/Q*-" (jlip . In the ID periodic 
setting there are only two nonzero components of Up. For atoms a and (3 the term £99 {{ua — up) /e) 
gives the following contribution to the right-hand side vector: 

(p'{{Ua - Up)/e) (Cfc - Up) , 
and the following contribution to the stiffness matrix: 



1 
e 

1 



ip"{{ua - up)/e) {ekel - ekU^ - U pel - UpU^) 



= -(p"{{ua - up)/e) (cfc - Up) {sk - Up)'^ . 

OCP 

Here is the unit vector corresponding to Uj^.. Adding contributions of £"3 to the previously 

OOP OOP 

computed contributions of Ej^ and E2 concludes assembling of the stiffness matrix and the 
right-hand side vector. 

The described implementation methodology of the QCP method is not confined to the periodic 
ID setting, but can easily be generalized to 2D or 3D, triangles whose sides are not aligned with 
the atomistic lattice, non-triangular mesh, etc. For that, one just needs to go through all pairs 
of atoms and sum their contributions in the following way. Consider two atoms a and /?, whose 
positions are represented as 

Ua = Ul [u^J^ and up = U^ [uij^ . 
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Then the contribution of the term e<p ((uq — up)/e) to the stiffness matrix and the right-hand side 
vector can be easily expressed through the difference of their representation vectors Ua — Up. One 
should take into account that the vectors Ua and Up are sparse and implement the respective 
vector and matrix operations efficiently. 

However, going through all the atoms in the summation is expensive and should be avoided. In 
order to optimize the summation procedure, the following observation is to be made. For all atoms 
a and /? within one element, all contributions to the right-hand side vector and the stiffness matrix 
can be computed once per element. It is commonly implemented by choosing a representative atom 
well inside the element and computing its energy contributior0. Similar optimization can be done if 
the element boundaries (in 2D or 3D) are aligned with the atomistic lattice. Then one can identify a 
fixed number (which depends on the interaction distance and the lattice structure) of representative 
atoms at element edges and faces. Through these representative atoms all contributions given by 
atoms in different elements, in principle, can be computed once per mesh edge or face. Thus, in 
such case the number of operations would only depend on the number of elements K and would 
not depend on the actual number of atoms N which can be very large. However, if the element 
boundaries in a 2D or a 3D lattice are not aligned with the mesh, then the number of operations of 
the QCP method will be higher than 0{K). In such cases one may consider employing some cluster 
summation rules which approximately sum all the contribution of atoms close to the boundary with 
0{K) operations. 

4.5 QCP Method and Ghost Force 

In this subsection we show that the QCP method does not exhibit any ghost force whatsoever. It 
should be noted that absence of the ghost force was first noticed by Rudd and Broughton in [25] 
for the so-called coarse-grained molecular dynamics (CGMD) method in a general setting. The 
QCP method is essentially a zero-temperature static CGMD method, therefore the argument of 
|25) can be applied to it, with the modification that the argument of t25j concerned the linearized 
model whereas we consider a general nonlinear model. Here we adopt the Rudd and Broughton's 
argument to the QCP method in the ID periodic setting as well as in the general setting. 

Consider the zero external force fi = 0. Then obviously a uniform lattice Ui = ezi satisfies the 
equilibrium equations ([2]), which means that the total force on any atom is zero. Let ifc (1 < ^ < K) 
define the nodal atoms for the QCP approximation. Obviously a uniform lattice Ui = ezi satisfies 
the QC constraints ([6]), therefore it is an admissible QC lattice. We say that a certain QC method 
does not exhibit the ghost force if Ui = ezi is indeed the QC solution to the problem. Otherwise 
we say that the QC approximation introduces some ghost forces that cause the uniform lattice to 
lose its equilibrium. The following theorem states that QCP does not exhibit the ghost force. 

Theorem 1 (absence of ghost force for the ID periodic setting). Consider the zero external force 
fi = and a uniform lattice «•* = ezi. Let ij. (1 < k < K) define the nodal atoms for the QCP 
approximation. Then it*^'*^^ = [u^ ]k is the QCP solution, i.e. 



u, 



,0 

■'fc 



(l<A:<i^). 




^In a number of works the representative atoms where chosen at the element vertices, which has been shown 
recently to yield inaccurate results, see [18] for further discussion. 
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Proof. First, notice that in the definition of IT (jl3p the term JJ^^u'-''^^ equals to u'-', since the 
linear interpolation operator [/Q^ on the uniform lattice is identical. Second, in the ID periodic 
setting the uniform lattice is in equilibrium: 



which concludes the proof of the theorem. 

We can observe that in the proof we essentially used only the fact that the uniform lattice is in 
equilibrium and at the same time is identity under the interpolation operator U^'^. This leads us 
to the following theorem in the general setting. 

Theorem 2 (absence of ghost force for the general setting). Let 11 = E'tot + -^ext be the potential 
energy of a general atomistic system described by its degrees of freedom u G M^. Let the lattice vP 
be uniform in the sense that it is in equilibrium w.r.t. H (in other words, is a critical point of 
H). Let the QC approximation be chosen so that there exists a QC vector u'^'^^ corresponding 
to the uniform lattice = [/QC.j^o,QC_ j^/jg^ ^o,QC ^ qq solution. 

The proof of this theorem involves straightforward use of theorem conditions to show that (|23p 
also holds in the general setting. 

In fact, a result somewhat stronger than Theorem [2] can easily be proven: if in addition is 
a local minimizer of 11, then u^'Q^ is a local minimizer of n*^^^. This relies on a simple fact that 
if a functional is restricted to a subspace which contains a local minimum u^, then uP remains a 
local minimum of the restricted functional. 

5 Numerical Results 

First, we present the results of ID tests in the same periodic setting under which we did analysis of 
the QCP method. Then we present the results of a 2D test. In all tests we chose the Lennard- Jones 
potential ip{z) = z~^'^ — 2z~^ with a cut-off radius of 3.25. Such potential involves third nearest 
neighbor interaction in one dimension. The results by the four methods, namely: the projection 
method (QCP), the classical quasicontinuum method also known as the energy-based quasicon- 
tinuum method (QCE), the quasinonlocal quasicontinuum method (QNL), and the geometrically 
consistent reconstruction-based method (GCR) are presented and compared. For the GCR method 
we took the particular instance of the method with coefficients given in table I of [9] for the ID 
tests and from table II for the 2D test. 

5.1 ID tests 

5.1.1 Test with Localized External Force 

We considered a periodic lattice with a period of length 1 with N = 10 000 atoms: 



dU 

dui 



= 0. 




= 



(23) 



Ui+N = Ui + 1 ( — oo < i < oo). 
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Figure 3: Numerical error for ID test with no bulk force. Error is computed in a discrete W^'°° 
norm, m is the number of nonlocal atoms near the defect. 



To model a defect of the atomistic material in ID, we considered the external force to be zero for 
all atoms except atoms N/2 and N/2 + 1: 

-1 i = N/2 
1 i = N/2 + 1 
otherwise. 

Correspondingly, the nonlocal region was chosen to be comprised of the following m atoms (only 
even values of m were chosen): 

= [N/2 - m/2 + 1, N/2 + m/2]. 

Since the external force is zero in the local region, we should expect the interatomic distances not 
to vary essentially in the whole local region and therefore no additional nodal atoms are required 
inside the local region. Thus, the nodal atoms were chosen as 

ii=N/2-m/2 + l, 12= N/2-m/2 + 2, i^ = N/2 + m/2. 

The dependence of the error in numerical solutions u^'~^^, u^^^, u*^*-"^, and u^^^ on the length 
of the nonlocal region m is presented in figure [3l The error was computed as 

/..numerical „, exact \ ("„, numerical ..exact \ 



error = \u~''^''^ - -u^^^^*Ui,oo = max 

l<i<7V 

It can be seen that the errors of QCE and QNL do not converge to zero as m grows, although the 
error of the QNL method is generally smaller. Also, it can be seen that the errors of QCP and 
GCR converge to zero and are very close to each other, as expected from the above analysis. The 
convergence rate is exponential in m. The reason for the exponential convergence is that for the 
potential with a finite cut-off distance any disturbance of the uniform lattice decays exponentially 
(see, for example [6j where a QC approximation is analyzed in the case of linearized problem with 
second nearest neighbor interaction). 
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Figure 4: Numerical error for ID test with bulk force. Error is computed in a discrete W^'°° norm. 
DoF is the number of degrees of freedom which equals to the total number of nodal atoms minus 
one. 



5.1.2 Test with Non-Localized External Force 

In the second test case we chose the same periodic lattice with = 10 000, but the external force 
was taken to be the sum of the irregular and the regular component: 

f = f" + r^ 

where 

10 i = N/2 
-10 i = iV/2 + l 
otherwise, 

and 

The irregular component models the defect of material and the regular part models the external 
bulk force. In this case the deformation in the local region is not constant and hence the nodal atoms 
should be chosen in the local region as well. In the present test the nodal atoms were comprised of 
a number of equidistantly spaced atoms in the local region and all atoms in the nonlocal region. 

The dependence of errors in numerical solutions u^'^^, u^^^, u^^^, and u^'^^ on the number 
of degrees of freedom (DoF) is presented in figure HI The number of degrees of freedom is the 
number of nodal atoms minus one (because we fix the position of one atom to avoid uncontrolled 
shift of the lattice as a whole). For each value of DoF and each method we chose the length m 
of the nonlocal region in such a way that it minimizes the error for the fixed DoF. It can be seen 
that the QCP method and the GCR method again give very close results. The errors of u^''^^ and 
.^GCR g^j.g ggg]2 to decay with the first order in DoF. The errors of QCE and QNL are seen to decay 
for small DoF and be approximately constant for large DoF because of dominating interface error. 
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Figure 5: 2D test with a point singularity, actual configuration 




Figure 6: 2D test with a point singularity, reference configuration 
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Figure 7: 2D test with a point singularity, a typical triangulation 



5.2 2D test 

We now present a numerical test not under the one-dimensional periodic setting but under a general 
setting of the atomistic problem with defects. We considered the configuration of atoms shown in 
fig. [5] with the reference configuration shown in fig. [6l Such configuration models a point defect in 
a 2D lattice. It corresponds to a dislocation in the respective 3D lattice (i.e. if one creates a 3D 
lattice by replicating these 2D layers, then the point defect will turn into a line defect which is an 
edge dislocation). The lattice was comprised of A'totai = 3365 atoms. The behavior of the atomistic 
solution resembles the point singularity: the strain decays as where r is the distance to the 
defect (the asymptotics of the strain is a general behavior of the strain for the edge dislocation 

m- 

The boundary conditions were set as follows. First, the full atomistic problem with A'totai = 3365 
atoms was solved using the Neumann boundary conditions (i.e. the positions of the atoms on the 
boundary were not fixed). This solution was used to fix the coordinates of the three outmost layers 



of atoms (shown as smaller circles in fig. 5(a) and 6(a)). Thus, the Dirichlet boundary conditions 
were used in the computations with three layers of atoms fixed. This was done in order to eliminate 
the effect of the boundary conditions on the lattice near the boundary and thus to eliminate the 
effect of the free boundary on the error of the QC solution. The typical triangulation used is shown 
in fig. [71 Notice that the three outmost layers of atoms do not participate in the triangulation. In 
fact, they also do not participate in the atomistic solution, only yielding the external force on the 
boundary atoms. Thus, the effective number of atoms in the computations was N = 2789. 

By similarity with the ID tests, the error of numerical solution was computed in a vector 
Ty^'°°-norm, which is defined in the following way. We loop over all triplets of neighboring atoms 
in the numerical and the exact solution and compute the Jacobian of the mapping of the triangle 
(corresponding to these three atoms) in the exact solution to the corresponding triangle in the 
numerical solution. The 2-norm of the difference of this Jacobian with the unit matrix was taken 
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Figure 8: Numerical error for the 2D test. Error is computed in an analog of W^''^ norm. The 
dashed vertical line corresponds to the total number of degrees of freedom in the atomistic system. 



as the error for one triangle. Then we take the maximum over all such triplets of neighboring 
atoms, which is the sought W^'°°-norm.. 

The dependence of the error in numerical solutions by different methods on number of degrees 
of freedom (DoF) is presented in figure El It can be clearly noticed that the errors of the QCP 
method are much less than those of the other methods and the error of the original QCE method 
is the highest. We argue that the reason for the error of the QCP method to be much less than the 
errors of other methods is that it handles interfaces between the mesh triangles in the local region 
in a more accurate way. To verify this, we also implemented the modified QCP method (QCPm) 
which is the same as QCP method everywhere except at the interfaces between the mesh triangles 
in the local region, where the standard QC discretization was used. In other words, QCPm uses 
the QC discretization everywhere in the local region except at the interface between the nonlocal 
and local region, where QCPm behaves as QCP. The results on fig. [8] show that the errors of the 
QCPm method and the GCR method are very close to each other. It can also be seen that errors 
for GCR and QCPm remain constant for small and moderate DoF and start to decay only when 
the DoF of the numerical solution is close to DoF of the original atomistic problem (the total DoF 
of the atomistic problem is shown with the dashed line in fig. [8|) . On the contrary, the error of the 
QCP method decays steadily as DoF increases. 

6 Discussion and Conclusion 

The presented QCP method coincides with the classical QC approximation for nearest neighbor 
interaction (n = 1). It is not equivalent to the QNL method, although both methods effectively 
eliminate the ghost force for second nearest neighbor interaction (n = 2). As compared with the 
GCR method, we may conclude that QCP is different from the particular instance of the GCR 
method presented in but there exist other instances of the GCR method which are equivalent 
to the QCP method (more precisely, to the QCPm method) in the way how these methods treat 
the local-nonlocal interface. 

The results of ID and 2D numerical tests show the following patterns. The error of the QCP 
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method was always lower than the errors of QCE and QNL. The error of the QNL method was 
lower than that of the QCE method, although both errors do not converge to zero. 

The error of the QCP method was essentially the same as the error of the GCR method for ID 
tests and lower for the 2D tests. We argue that the reason for such behavior is that in ID case 
the disturbance of the uniform lattice decays exponentially, whereas the disturbance in the 2D test 
decays much slower: as the inverse of the distance to the defect. This causes large interface errors 
between elements within the local region for all methods except QCP since the former use the 
Cauchy-Born extrapolation. Not using the Cauchy-Born extrapolation in the local region results 
in a somewhat larger number of operations for assembling the stiffness matrix and the right-hand 
side vector, but may give much more accurate results as compared to the other methods. In case if 
such increased accuracy is not essential then one can use the modified QCP (i.e. QCPm) method 
which treats the interface between nonlocal and local region as the QCP method and uses the 
Cauchy-Born extrapolation in the local region. The QCPm method yields essentially the same 
results as the GCR method but, as we believe, is easier to implement, analyze, and generalize to 
other problems. 

GCR is a cleverly designed method and gives a very general criterion to remove the ghost force at 
the local-nonlocal interface. However, based on the results of analysis and numerical experiments we 
infer that the QCP method (including QCPm, which is a simpler but less accurate modification of 
the projection method) is easier to implement and analyze, does not depend on any free parameters 
(such as reconstruction coefficients), and may be even more accurate than GCR in general, in spite 
of the fact that QCP may be expressed as a particular case of GCR (with a difference that they 
treat the interface between elements in local regions differently). First, as discussed in subsection 
14.41 for implementation of QCP one just needs the atom position's representations through the 
nodal atoms. No other geometric information is needed for the implementation, which makes 
it easier than GCR, since the later requires the a priori tabulated coefficients of reconstruction 
and needs to determine distance to the local-nonlocal interface for each atom near the interface. 
Second, the QCP method is also easy to analyze: one can benefit from the well-developed theory 
of finite elements which offers a powerful method of reducing the problem of convergence to the 
problem of approximation (see e.g. [16\ I22j). Third, the QCP method does not require determining 
any problem-dependent parameters beforehand. It is therefore more flexible: one needs to do less 
investment for solving the new problems. For example, if one wants to apply the QC method for 
metallic alloys, whose atomistic lattices are not uniform due to presence of atoms of different metals 
in the lattice, one just needs to describe this lattice in terms of reconstruction of the non-nodal 
atoms through the nodal atoms. It should also be noted that the QCP method is essentially a 
particular case of the CGMD method originally proposed for coarse-graining the finite temperature 
molecular dynamics [231 [Ml [25] . Therefore it may potentially have wider applications. Fourth, as 
we have seen from the 2D numerical test (subsection 15. 2p . the QCP method may be more accurate 
since it also offers a natural way to eliminate the edge error in the nonlocal region as compared to 
the standard Cauchy-Born extrapolation. 
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